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Abstract. An existing phase-field model of two immiscible fluids with a single 
soluble surfactant present is discussed in detail. We analyze the well-posedness 
of the model and show that it is mathematically ill-posed for a large set of 
physically relevant parameters. As a consequence, critical modifications to 
the model are suggested that substantially increase the domain of validity. 
Carefully designed numerical simulations offer informative demonstrations as 
to the sharpness of our theoretical results and the qualities of the physical 
model. A fully coupled hydrodynamic test-case demonstrates the potential to 
capture also non-trivial effects on the overall flow. 



1. Introduction 

The presence of surface active substances may greatly affect the physical prop- 
erties of fluid mixtures. Indeed, these effects are used critically in many important 
applications in everyday life; detergents and oil-water emulsions in food are two im- 
mediate examples. The fact that surfactants lower the surface tension is exploited 
in both of these cases: detergents make the water more "wet" , and an emulsifying 
agent stabilizes an emulsion by preventing small droplets to coalesce. 

For humans, probably the most critical everyday usage of surfactants is made in 
the alveoli in the lung, where pulmonary surfactant, amongst other things, prevents 
lung collapse at the end of expiration [22, 26]. 

An interesting and very striking realization in vivo was reported recently in [32] 
where chemotaxis was implemented for small droplets of fluid in a bulk solution, 
physically contained in a maze. The net transport at the millimeter-scale, and 
the subsequent solution to the maze, was achieved by a clever usage of surfactant 
and a pre-existing pH-gradicnt. Quite likely, such a constructive set-up could find 
applications in lab-on-a-chip manufacturing. 

Consisting of hydrophobic "heads" and hydrophilic "tails" , surfactant molecules 
have a strong preference to occupy sites at the water-fluid or water-gas interfaces. 
Below the critical micelle concentration (CMC), surfactants therefore adsorb ef- 
ficiently to the interfaces where their physical effects become prominent. Above 
the CMC, additionally, spontaneous formation of stable groups of surfactants - 
micelles — occurs in the bulk solution [38]. 
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Given the 'thermodynamical' signature of surfactants; that is, the diffusion- 
limited flow into the interfaces and the spontaneous creation of highly regular mi- 
celles from an unordered state in the bulk, modeling through some kind of system's 
energy assumption is a tempting approach. 

Ariel, Diamant, and Andelman [1, 9] have successfully postulated free-energy 
terms with theoretically convincing properties. Their approach is inherently sharp 
in that they set up equations for the interface, the sub-surface region, and the bulk. 
A later work [38] shows that this methodology can be extended to include also the 
region above the CMC. A great feature with this type of modeling is the fact that 
all properties of the model result from a single postulated entity; the system's free 
energy. 

Given the multitude of scales present and the complex coupling to hydrodynam- 
ics, numerical experiments become important as tools to gain a better understand- 
ing and a fuller physical insight. 

In a so-called sharp interface method, the interface is considered to be infinites- 
imally thin, and its exact location is represented either explicitly (front-tracking 
with e.g. Lagrangian markers), or implicitly (e.g. level-set [55] and volume of fluid 
(VOF) methods [28]). The evolution of the surfactant concentration on each in- 
terface can be described by a partial differential equation on this time-dependent 
manifold [47] . Techniques for solving this PDE and include insoluble surfactants in 
multiphase flow simulations have been developed based on several different interface 
representation techniques [3, 28, 31, 33, 42, 55]. 

When the surfactants are soluble also in the bulk, source terms due to adsorption 
and desorption terms enter this PDE, and this must be coupled to a PDE for 
the bulk concentration of surfactants with appropriate boundary conditions for 
the surfactant flux. See e.g. Eggleton and Stebe [13] for an early reference. To 
simplify matters for simulations, it is in this reference assumed that the adsorption 
and desorption from the bulk is diffusion dominated and that the bulk surfactant 
concentration is spatially constant. Work to consider the full problem has started 
only recently and is currently an active area of research [5, 30, 39, 50, 57]. 

Given a thermodynamic description as a postulated free energy of the system, 
including surface effects introduced via gradient energy terms [34, 35, 36, 44, 51, 52], 
it is attractive to represent the interface as a rapid but continuous transition in a 
concentration variable. This is the essence of phase-field or diffuse interface model- 
ing. Without the need to explicitly track the interface, even complicated topological 
transformations can be studied, and the coupling to the full hydrodynamic flow is 
fairly straightforward [27]. 

However, this has only recently begun to be explored for surfactant laden inter- 
faces. In [44], and with some additions, in [36], a seemingly natural diffuse interface 
model for surfactant laden liquid-liquid mixtures is presented, and solved with lat- 
tice Boltzmann methods. The model is verified against classical results, and some 
cases of convection dominated flow, such as a deforming droplet in a shear flow 
are presented. This formulation is promising from a computational and physical 
point of view, but it is rather complex and its mathematical properties have not 
been investigated. As we will show below, this model is in fact mathematically 
ill-posed under certain conditions. The purpose of the present paper is to clarify 
the conditions under which the model is ill-posed, and also to present variants that 
behave better. 
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In Section 3 we thus present the analysis of the model, and we also derive three al- 
ternatives in the form of modifications to the original model. The properties of these 
alternatives are investigated in Section 4 using a highly accurate one-dimensional 
numerical scheme, and in a more qualitative sense in two spatial dimensions in 
Section 5. An outlook with conclusions is finally found in Section 6. 



2. Diffuse interface model for surfactant flow 

In this section we shall discuss our 'baseline' phase- field model for multiphase flow 
incorporating surfactants; — later we shall have reasons to alter the actual model 
by incorporating new terms. Our starting point is a model originally presented in 
[44], but incorporating the work also of others. 

Specifically, a diffuse interface model for surfactant-controlled emulsification was 
proposed in [52] and a related model for oil-water-detergent mixtures appeared in 
[34] . A model with more easily understood interface adsorption properties was thus 
proposed in [44], incorporating free energy terms developed previously for a sharp 
model [10], but drawing also on some earlier work [35]. For clarity we take a non- 
dimensional version as our baseline model (see Appendix A for the precise relation 
to the original model in [44]). 

The degrees of freedom are the phase-field variable <p G [—1,1], the surfactant 
volume fraction ip € [0, 1], and the fluid velocity field u. The governing equations 
for cj) and tp take the form of Cahn-Hilliard-type equations [27], 

(2.1) ^ + V . W = ^-V.^V W = 1% 

dtb 1 

(2.2) + V ■ (V>u) = — V-AfyV/v, 

in terms of Peclet numbers Pe, mobilities M, and chemical potentials p which will 
be prescribed below. As indicated in (2.1) the model is somewhat simplified in that 
it assumes the "shallow quench" limit [41, Sect. 9], and hence that the mobility 
for tf> is assumed to be constant (and in fact normalized to unity by scaling 
the Peclet number Pe^ appropriately). This simplification is quite standard and is 
adopted here mainly for convenience. 

Eqs. (2.1)-(2.2) are coupled to the Navier-Stokes equations in the form 

(2.3) V-u = 0, 

(2.4) p(^+u-Vu)= -VP + ^V ■ (MVu + (Vu) T ]) 

1 



CaCnRc 



(0V/X0 + ipVp,^) 



where the pressure tensor P enforces the incompressibility condition (2.3), and 
where Re is the Reynolds number, Ca the Capillary number, and where, respec- 
tively, p and v denote the fluid's density and kinematic viscosity. In the general 
case, p and v may depend on <j) as well as on other state variables, requiring that 
additional effects be taken into account [11]. For our present purposes and for 
simplicity they are assumed to be constants. 
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2.1. The Ginzburg-Landau free energy. Recall that in diffuse models of binary 
fluids in equilibrium, the phase-field variable <j) typically describes a planar interface 
at the origin of the form 

(2.5) <j>(x) =tanh(x/Cn), 

where Cn is the Cahn number expressing the ratio between the interface width and 
the characteristic length scale. 

In (2.1)— (2.2) the chemical potentials are derived from a Ginzburg-Landau free 
energy functional, 

(2.6) f Fdx= [ F^+F i> +F 1 +F UK dx. 
Jo. Jn 

The classical constant mobility Cahn-Hilliard potential [41, Sect. 4] is given by 

(2.7) F , = -f4^(Vtf, 

where we note that the polynomial part is a "double well" potential, [(</> 2 — l) 2 — 1]/4, 
and hence expresses a preference to pure (non- mixed) phases <f> = ±1. Also, the 
explicit penalty term for steep gradients causes any interfaces to become diffused. 

For "0, a logarithmic free energy is rather preferred as this provides for certain 
analytical properties such as possessing an isotherm relation as outlined in Sec- 
tion 4.2, 

(2.8) Fy = Pi [i/> log i/> + (1 - V) log(l - if>)] . 

The temperature-dependent constant Pi (cf. (A. 11)) takes the role of a diffusion 
coefficient for tjj and Fy is thus the energy potential governing the entropy decrease 
of mixing the surfactant with the bulk phase. We shall later find reasons to modify 

(2.8) since it, unlike (2.7) does not contain a square gradient term. Note that a 
logarithmic free energy is to be combined with a degenerate mobility which vanishes 
at the extreme points ip S {0, 1} [41, Sect. 4.2]. The usual choice is simply My — 

— ip) which will also be used here. A point in favor of formulating the Cahn- 
Hilliard equation using this kind of mobility is that, unlike the constant mobility 
case, it can rigorously be shown to produce solutions < tp < 1 [41, Sect. 4.2]. 

Surfactant molecules tend to move to a point and orient themselves in such a way 
that they are at ease. Owing to their nature this generally means at the interface 
between fluids, with the hydrophilic part pointing towards a fluid of the system 
and the hydrophobic part pointing away from this fluid. In the current model, F\ 
is the surface energy potential accounting for this adsorption, 

(2.9) F X = -^(V0) 2 . 

This particular form is inspired by the corresponding term in the sharp interface 
model of [9, 10]. The square gradient acts as a diffuse version of the sharp inter- 
face indicator function and can rigorously be interpreted as a nascent Dirac delta 
function as outlined in Section 3.2.2 below. 

Finally, to conclude the specification of the free energy, the term _F CX takes the 
form 

(2.10) F - = -m*? 
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and penalizes free surfactant in the respective phases. Similarly to F^, f ex is an 
enthalpic term measuring the cost of free surfactant. However, due to the factor 
<fi 2 , F cx is inactive at an interface where ^0. Actually, (2.10) is just the simplest 
form among many possibilities (cf. [34, 35, 52]). To some extent, (2.9) and (2.10) 
are complementary: F\ locally attracts surfactant to an existing interface while F ex 
globally counteracts the occurrence of free surfactant. Finally, and as we shall see 
in Section 4.1, (2.10) can also be regarded as defining the bulk solubility. 

Through variational derivatives of the free energy F with respect to the two 
variables <p and ip one now obtains the chemical potentials 

(2.11) = ^ = -4> + (f> 3 - *=^A0 + ^VA0 + ^ttW • + 2^M. 

,,12, W = ^=Pi.og I ^- C ^(V^ + ^. 

Note that, carrying out the differentiation in the right-hand side of (2.2) for the 
logarithmic term in (2.12), we get 

^ F^ V -^ VPi ( l0g l^)-F^ A ^ 

provided that < ip < 1. 

This completes the specification of the gradient flow (2.1)-(2.2) which together 
with the Navier-Stokes equations (2.3)-(2.4) and suitable boundary conditions make 
up our basic multiphase-surfactant model. In what follows we shall initially be 
concerned mainly with the diffusion-controlled part (2.1)— (2.2) (with u = 0) and 
postpone computational experiments with the full hydrodynamic set of equations 
until Section 5. Testing (2.1) with and (2.2) with /Lt^,, respectively, we find under 
natural boundary conditions that 

(2.14) l/^^ = -/ fi p^l|VM 2 + ^l|V^|| 2 ^<0 

for > 0. This relation expresses the decay of the free energy and is required for 
thcrmodynamical consistency. For the full hydrodynamic set of equations we find 
similarly by testing (2.4) with CaCnReu the general energy identity 

d f pCaCnRe, _ , 

' ' \u\ 2 +Fdx = 



dtJ Q 2 

(2-15) - / -L||V^|| 2 + ^||V^|| 2 + ^^||Vu+(Vun| 2 ^. 

3. Analysis 

In this section we first show that the free-space PDE defined by (2.1)-(2.2) with 
u = is ill-posed in relevant parameter ranges. We then return to the baseline 
model (2.1)-(2.4) and modify it in order to provide for well-posedness. Three 
different modifications are suggested, all with differing motivations and properties. 

3.1. Ill-posedness. We thus consider the case u = in (2.1)-(2.2) for a one- 
dimensional free-space formulation. We make the quasi steady-state ansatz <f> — 
<^oo + Su and ip = ipoo + Sv, where (poo/ipoo are equilibrium solutions as t — > oo 
and where 8 is a small parameter. Assuming the existence of equilibrium solutions 
is motivated by the fact that (2.1)-(2.2) describe a gradient flow and naturally 
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possesses (non-unique) equilibria [48, Sect. 4]. — Fundamentally, the lack of any 
equilibrium would violate the second law of thermodynamics. 

Upon linearization in S and keeping only the principal part of the operator we 
get after some work the linear system of PDEs 

-^oo n 4 Cn 2 £>(0oo) 

^ 2 
)D{<t>o 



(3.1) 







v t 





Cn" 



_Cnl 
2 

^oo(l- 



D 



Pe^ 

Pe 





u 




V 



where D = d/dx. Without further ado we take the "frozen coefficient" Fourier 
transform of (3.1) by simply exchanging D — > —iuj, 



(3.2) 







v t 





Cn' 



2 

V-=o(l 



Pej, 

-y. oo )n(0 oc 



Cir 
2 



Pc tf 



Pi 

_^p 

Pe 



2 



Call the above matrix A. Firstly, since the characteristic equation is real, complex 
eigenvalues must come in conjugate pairs. Secondly, from a quick calculation we see 
that det A cx Pi — Cn 2 /2-ijj oc> (D((j) rx> )) 2 and hence, when the determinant is negative 
there is a positive real eigenvalue indicating some kind of instability. Slightly more 
elaborate calculations are now required to show that there is in fact an eigenvalue 



(3.3) 



A, 



1 

Pe^ 



Cn' 

~2~ 



<M£(0oc)) 2 



Pi 



•0(1), 



while Ai is asymptotically negative and scales as O (w 4 ) . From (3.3) we see that 
A2 may well be positive for w sufficiently large. If this is the case the PDE violates 
the Petrovskii condition which is necessary for well-posedness [19, Theorem 4.5.2]. 

Using results to be developed in Section 4 we now derive some simplified con- 
ditions for this unphysical instability to occur. Let us consider it inherent in the 
model that both ip^ and (-D(^oo)) 2 are largest at the interface (taken to be at 
x = 0). Then ipoo < "0oo(O) = ipo, the surfactant loading at the interface and, 
provided that the interface sharpness is independent of the surfactant loading, 
(£)(0oo)) 2 < (£>(</>oo(0))) 2 ~ 1/Cn 2 (see Section 4.1). Using the isotherm rela- 
tion ipo r~j ipt,/(ip c 4- ip b ) (again, see Section 4.2 below), an a priori approximate 
sufficient condition for instability is therefore 



(3.4) 



Pi < 



expressed in terms of the bulk value tpi, — ^^(±00) and the Langmuir adsorption 
constant ip c discussed thoroughly below. Using the relation (4.15) for the latter 
this is equivalent to 



(3.5) 



2 Pi 



1 - 2 Pi 



X 1p c = 



2 Pi 



1 - 2 Pi 



x exp 



1 + 1/Ex 
4Pi~ 



provided that Pi < 1/2. That is, sufficiently large values of ipb lead to an instability 
near the equilibrium solution (which therefore does not exist). 

This situation is somewhat remindful of results for a close relative of the Cahn- 
Hilliard equation, namely the Thin film equation. This model can also be derived 
from a Ginzburg-Landau energy which, as in (2.14), by construction decays with 
time. Through a clever constructive argument, however, one can show that for 
a sufficiently large total mass, there is a smooth initial profile for which finite- 
time blowup in the form of indefinite focusing occurs (see [54] and the references 
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therein). Unfortunately, it seems difficult to transfer these types of arguments into 
the current setting. 

Example. In several numerical experiments in [44] using the lattice Boltzmann 
method, Pi € [0.08, 0.20] and Ex ~ 1. To be concrete and still use actual (published) 
values, we take Pi = 0.1227 with Ex = 1. Then (3.5) becomes tp b > 5.526 x 10~ 3 
and under this condition we expect the model to be ill-posed. Since many of the 
experiments in [44] (see for example Fig. 1 and 2) are made in this regime, these 
results are misleading at the very least. One cannot help speculating that the lattice 
Boltzmann method used enjoy "hidden" higher order derivative terms that help 
increasing the stability. In Figure 3.1, p. 10 these predictions are illustrated using 
the numerical method developed in Appendix B. We note that similar experiments 
are made anew in [36] under only slightly altered conditions. Here, the parameters 
in the simulations appear to have been kept within the limits for well-posedness. 
Using this model substantially limits the physical parameter space that can be 
explored. 

Incidentally, (3.5) disqualifies the label "numerical stabilization term" for the 
bulk-energy F cx in (2.10) as used in [44, p. 4] and again in [36, p. 9167]. Increasing 
their variable W in (A. 11) (that is, "increasing" the numerical stabilization) is 
the same thing as decreasing Ex, thereby driving more surfactant to the interface. 
Assuming that the other parameters Cn and Pi are kept fixed, according to (3.5) 
this means that less amount of surfactant implies mathematical (as opposed to 
numerical) instability. 

Finally, and from the perspective of physical modeling, one could in principle 
accept the ill-posedness as an artifact of the fact that we are dealing with an 
artificially diffuse interface. However, in such a case it is natural to demand at the 
very least that the limiting problem Cn —> is stable. For the eigenvalue A2 in 
(3.3) to be negative as Cn — ► one clearly has to give up the natural requirement 
that the interface sharpness is independent of the surfactant loading. Indeed, our 
approximate bounds (3.4)— (3.5) are fully independent of Cn as a reflection of this 
assumption. 

3.2. Three alternatives. In this section we seek viable alternatives to the model 
just analyzed. Our first modified formulation involves using the full form of the 
logarithmic free energy in (2.8). This form contains a square gradient term and 
leads to an additional fourth order dissipative term in the PDE for the surfactant 
concentration if). The downside is that the resulting model has more complicated 
mathematical properties and we therefore look for simpler but still reasonable al- 
ternatives. Two such suggestions are obtained by using a derivative-free diffuse 
interface Dirac delta function when defining F± in (2.9). 

3.2.1. Completing the logarithmic free energy. Rather than (2.8), the complete log- 
arithmic free energy is actually given by ([41, Sect. 4.2, Eqs. (20), (21)]; see also 

[14]) 

(3.6) = Pi [01og0 + (1 - V) log(l - 0)] + - VO + ^p(W) 2 , 

with a a new parameter. Accordingly, our "Model 1 " is obtained by exchanging 
the definition of F^ in (2.8) with (3.6). While the chemical potential fi^ in (2.11) is 
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not affected by this addition to the free energy, fi^ in (2.12) is to be replaced with 

(3.7) H - P.log^ - <^ (Vtf - 9f AO - ^ + ^ 

Evidently, one more boundary condition than before is needed since now the equa- 
tion for tp is fourth order. 

A possible explanation to the fact that the square gradient term is missing in 
the original presentation of the model in [44] is that it is also omitted in the earlier 
references [9, 10, 35]. In [9, 10] we note that the term has no meaning since the 
interface model there is sharp. Interestingly, the square gradient term is in fact 
included in the surfactant models found in [34, 51, 52], but then for the constant 
mobility case. As remarked in [44] (and as will be evident in Section 4.2 below), 
such models do not possess a realistic adsorption isotherm. 

This also seems to be the place to mention how the two order parameters (<p, tp) 
are to be interpreted. Rather than a true multicomponent system as treated in [14], 
we are dealing with small, often very small, concentrations of surfactant. A kind 
of two-step approximation is therefore employed where a Cahn-Hilliard equation 
with degenerate mobility governs the surfactant (tp) and "non-surfactant" (1 — ip) 
phases. In turn, the non-surfactant phase (e.g. oil and water) is treated by a classical 
constant mobility Cahn-Hilliard equation in the phase-field variable (p. The three 
concentration components of the flow should therefore rightly be understood as 

(3.8) [il>, (1 - ^0(1 + 4>)/2, (1 " ~ 0)/2] 
and are thus approximated with 

(3.9) [^,(l + 0)/2,(l-0)/2] 

so that the volume is only preserved up to O (tp) . 

We conclude by offering a few comments on the lateral interaction term —a/4-tp 2 
in (3.6). This term is included also in [9, 10, 36] where it is shown to lead to the 
Frumkin isotherm (see Section 4.2). Interestingly, the same term appears also in 
[34, 35, 51, 52], but here with the opposite sign. Clearly, a positive value of a favors 
clustering of surfactant since the term can be interpreted as a free energy decrease 
for pairs of surfactant molecules. By differentiating twice the non-gradient part of 
(3.6) with respect to tp we further find that 

(3.10) ^ [Pi[V»logV+(l-^)log(l-^] + ^(l-^)] =— -|. 

Hence, for < tp < 1, this part of the free energy remains convex in tp provided 
that we choose a < 8 Pi. 

3.2.2. Gradient-free Dirac delta functions. We have already commented that the 
term - Cn 2 /4 • tp(Vcp) 2 in (2.9) results fr om using the square gradient as a diffuse 
version of the sharp interface indicator function. Indeed, for the planar interface 
equilibrium solution (p(x) = tanh(a;/ Cn), we have that (V</>) 2 = Cn -2 sech 4 (a;/ Cn). 
Over the real line, the quartic hyperbolic secant is a nascent Delta function in the 
sense that 3/(4 Cn) sech 4 (a;/ Cn) — > S(x) (convergence in distribution in the sharp 
interface limit Cn — > 0). 

Noting that the function 1 — tanh 2 (x / Cn) = sech 2 (x / Cn) also defines a nascent 
Delta function by virtue of the limit 1/(2 Cn) sech 2 (a;/ Cn) — > S(x), a tempting 
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replacement for (V</>) 2 in (2.9) is 2/(3 Cn ) • (1 — (f> 2 ) after appropriate scaling. 
However, as we shall see, this choice has a different behavior near the origin and 
it shall later be convenient to use the slightly altered scaling 1/ Cn 2 -(1 — <fi 2 ) (see 
Section 4.2). We thus define our "Model 2" by replacing (2.9) with 

(3.11) F 1 = -^(l-0 2 )- 

Taking variational derivatives we obtain the new chemical potentials 

(3.12) ^ = _0 + ^3_^_ A + 

(3.13) />,, Pilog ' ' - 1 



, 2 



1-ip 4 T 4Ex 



where it is clear that Model 2 can be implemented by simply skipping the term F\ 
altogether and substitute Ex — > 1/(1/Ex+1) in (2.10). 

Although the square hyperbolic secant is a well-known nascent Delta function 
(it is the derivative of the Fermi-Dirac function), in the present context one can 
actually continue to use the quartic hyperbolic secant, without introducing any 
derivatives. A simple replacement for the square gradient achieving just this is 
1/Cn 2 -(1 - (f> 2 ) 2 . Accordingly, our "Model 3" uses, in place of (2.9), 

(3.14) Fl = _I V ,(i_^ )2) 
with the associated chemical potentials 

(3.15) H = -4> + cf> 3 - + (1 - 4> 2 )U + 

,, 16) H . m ^-Q=£f^. X 

Evidently, Model 3 is stiffer than Model 2 and it will also be shown to produce a 
sharper equilibrium profile. Figure 3.1 and 3.2 below show representative sample 
simulations of all four models thus far considered. These results were obtained with 
the numerical method discussed in Appendix B. 



3.3. Sample simulation: ill-posedness. In Figure 3.1 we compare an unstable 
and a stable case of the original surfactant phase-field model (referred to as "Model 
0" from now on). The initial data, in the form of a uniform surfactant profile, 
was chosen in accordance with the stability criterion (3.5) and in such a way that 
both cases are in close proximity to the boundary of the region of well-posedness. 
During the simulation time displayed here it holds for the numerical solution tp 
that < i/j < 1 and that the associated Ginzburg-Landau (2.6) energy is decreasing. 
Shortly after the displayed simulation time the numerical solution becomes negative 
such that the free energy formally becomes multivalued. 

In Section 4.2 we more carefully evaluate the sharpness of the condition (3.5) 
after developing some more concepts (see Figure 4.1). 

Examples of all alternative models proposed in Section 3.2 are displayed in Fig- 
ure 3.2 were we simulate the ill-posed case from Figure 3.1 anew. Unlike Model 0, 
for this choice of parameters, all three new models are perfectly stable. 
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Figure 3.1. Time snapshots of the surfactant concentration for 
an unstable (solid) and a stable (dashed) case. All parameters 
are the same in the two cases (Cn = 1/6, Ex = 1, and Pi = 
0.1227), but the amount of surfactant differs by a factor of two. 
The initial data for ip(x) is a flat profile with height 0.012 and 
0.006, respectively, and tfi(x) is initially and in both cases set to 
tanh(ir/Cn). In the top two graphs the instability has not yet 
developed. In the two graphs below the O (w 2 )-part of (3.3) is 
positive inside the indicated small region near the origin, and the 
instability immediately becomes manifest. In the bottom graph 
the free energy (2.6)-(2.10) is plotted as a function of time t and 
the times for the four snapshots are indicated by circles. See text 
for further comments. 



4. Analysis and numerical experiments in ID 

In order to evaluate and compare the proposed models we now proceed to derive 
some analytical estimates, most of which we test through numerical simulations. We 
discuss equilibrium solutions, adsorption isotherm relations, and we also evaluate 
experimentally the diffusion-controlled adsorption dynamics at the interface. 

4.1. Planar equilibrium solution at constant surfactant concentration. If 

tp is held constant at a bulk value i/fy, steady-state of <j) implies for both our starting 
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model (2.11) and for Model 1 in Section 3.2.1 that 



(4.1) 



On' 



Acj)- 



Cn 



1 



= o. 



2 r 2 r " r 2 Ex 
Generally, steady-state requires /u^ = constant, but it is not difficult to see that a 
solution which is odd about the origin (say, <j>(±oo) = ±(1 — e)) must in fact have 
fj,$ = 0. Rewriting we get 



1 



1 



2 Ex 



(1_^)^-A0 



(4.2) 



V 1 ' 



ACn 2 
! 2 



A</> b "\ 



:Cn2 /2 



This is just a new Cahn-Hilliard steady-state equation in the scaled variable ip := 



b^cj) with Cahn- number Cn^, = Cn \/T~ 



From the classical solution ip — 



tanh(x/Cn v ) we therefore get 



(4.3) 



with 



4>(x) = ebb tanh I 4>b 



(4.4) 



= 1 - 



2 Ex 



At the interface we have for this solution that Cnd/dx <fi(0) = 1 — (1/ Ex— 1)^/2 + 
O (ipl) implying that the sharpness of the interface is independent of the surfactant 
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loading up to O (ipb) wnen Ex = 1. The bulk behavior as x — > ±00 is ±(pb which 
means that the coefficient Ex controls the bulk solubility. 

Proceeding in a similar fashion with Model 2 using (3.12) we obtain the equilib- 
rium solution 

(4.5) <t>{x) = 6 tanh (<fr>^) > 

i \ h 

Ex/ 2 ' 



(4.6) 



= 1-1 



and where Cn d/dx 0(0) = 4>i = 1 — O (ipb) and is < 1 for any positive value of Ex. 
Finally, for Model 3 using (3.15) we get after similar manipulations 

(4.7) 



4>(x) = <pb tanh I 4>b 



(4.8) 



1-1 



1 

2 Ex 



Cn^/T 

^ 06 



0b 



Not only are (4.7) and (4.3) identical in form, but (4.8) and (4.4) in fact also 
agree up to order O (tpf)- Consequently we have again that Cud/dx cp(Q) = 1 — 
(1/Ex— 1)06/2 + O so that for Ex = 1, the sharpness of the interface is 
independent of the surfactant loading up to O (ipf)- 

Conclusion. For all four models, the parameter Ex controls the bulk solubility. 
Model 0, 1, and 3 all agree closely in terms of the bulk value cpb and also in the 
sharpness of the phase-field interface <p = 0. By comparison, Model 2 has a more 
diffuse interface. 

4.2. Adsorption isotherm. We now turn our attention to equilibrium profiles for 
0, following the line of reasoning in [44] closely. Since at steady-state one must have 
that the chemical potential is constant, the basic approach is to solve the equation 
HMx) — fJ>tfi b f° r ^ where ipb as before is the bulk concentration 0(oo). We write 
the chemical potentials for all four models in the forms 

06 - 1 



(4.9) 
(4.10) 



Pi log 
Pi log 







B h 



B 



4 Ex 
1 



1-0 4Ex^ 

Subtracting and introducing the intermediate variable ip c { x ) we get t ne relation 

1 , ,0 



(4.11) 



Pi log ip c {x) = B - B b 



4 Ex 



2 ) 



in terms of which the steady-state profile is given by 

06 0b 



(4.12) 



ip{x) 



ipb + c (aO(l - 06) 
The difference B — Bb is given by, respectively, 



06 + 1pc(x) 



(4.13) 



B-Bh 



(V0) 2 

W) 2 



Cir 



2 ) 



A0- f(0-06) 



-i[(^_^)(2-^_ 02)] 



O (06) 



(Model 0) 
(Model 1) 
(Model 2) 
(Model 3) 
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Clearly, for Model 1, this line of reasoning needs to be augmented with additional 
assumptions or estimates in order not to be circular. For the other cases, with 
a fixed phase-field profile 4>{x) 4n, tanh(cc/ Cn) (as determined in Section 4.1), 
(4.12) yields a quite decent approximation as we shall see. 
Specializing x in (4.12) to the origin we get 

(4.14) ^ = _!£*_ 

Wb + Wc 

(4.15) mogj,c = -\(l+^)+Ofo). 

To arrive at (4.15), for Model we have to assume that Cn 2 (V(H0)) 2 = 1 + (V&) 
(cf. Section 4.1), while for Model 2 and 3 we only need to use the fact that </>(0) = 0. 

Eq. (4.14) is the Langmuir isotherm and ip c as defined by (4.15) is the Langmuir 
(equilibrium) adsorption constant. The fact that all models except Model 1 pos- 
sesses the same Langmuir isotherm makes them comparable and also explains our 
special choice of scaling when constructing Model 2 in Section 3.2.2. 

For completeness, we note that (4.12) for Model 1 becomes as x — > 0, 

(4-16) = . ~ , ^\ -ys+O^f,), 

tpb + ipc exp(-aVo)-R 

with i/j c still defined by (4.15), and where a = cr/(2Pi), R = cxp(- Cn 2 A^ /(2 Pi)). 
With R = 1, (4.16) is the Frumkin isotherm which essentially is an effect of the 
presence of the lateral interaction term — cr/4 • ip 2 in (3.6). Unfortunately, there 
is no evident relation between ipo and Aipo which can be used to close this line 
of reasoning. Therefore, for Model 1 the adsorption isotherm cannot be explicitly 
determined. 



4.2.1. Numerical isotherm. Using our one-dimensional spectral Galerkin code as 
outlined in Appendix B we have performed several numerical experiments with the 
adsorption isotherm relation. For all experiments in this section we used the same 
set of parameters: Cn = 1/6, Ex = 1, ip c G {0.0020,0.0056, 0.016,0.035,0.075} 
(with Pi determined from the relation (4.15)), and with ipt, sampled in the interval 
[10 _3 ,10 -1 ]. All parameters have been chosen to agree with those in [44, Fig. 2]. 
The simulations were started with the profile 4>{x) — tanh(x/ Cn) with ip(x) defined 
by (4.12) and, for convenience, <f>b = 1. Finally, for Model 1 we used the value 
<T = 8 Pi, obtained from the requirement that the non-gradient part of be 
convex (see (3.10)) and otherwise simply by trial and error to approximately match 
the visual appearances of the profiles for Model 2 and 3 (cf. Figure 4.4). 

In Figure 4.1 we numerically test the sharpness of our analysis of Model in 
Section 3.1. Clearly, the model is ill-posed for most of the parameters tested here, 
rendering the model very questionable. From the figure, it is also seen that the 
sufficient condition for ill-posedness (3.5) is quite sharp. 

In Figure 4.1 a numerical isotherm is obtained for Model 1. The Langmuir 
isotherm is not a bad model for small values of tp c , but the measured adsorption 
breaks off for larger values and also for higher surfactant concentrations. A few cases 
were found to be numerically unstable by producing an extremely sharp profile for 
</>. We have not been able to analyze this and do not know at present if this is an 
actual property of the model or a numerical artifact of some kind. 
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Figure 4.1. Attempt to verify the Langmuir isotherm (4.14) 
through numerical simulations of Model 0. This figure is intended 
to be an exact reproduction of the results in [44, Fig. 2]. How- 
ever, here we clearly see that the model is ill-posed as predicted 
by the theory in Section 3.1. Solid: the isotherm (4.14) for differ- 
ent values of ip c , circles: numerical values. The dotted line is the 
sufficient condition for ill-posedness (3.5) and crosses are used to 
denote (missing) unphysical solutions (e.g. large negative values). 




Figure 4.2. Results from simulations with parameters as in Fig- 
ure 4.1 but using Model 1 instead. Since there is no known adsorp- 
tion relation in this case, the Langmuir isotherm (4.14) in solid is 
displayed for reference only. A total of four unstable cases were 
detected (crosses). 
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1Cf 3 1Cf 2 1Cf 1 10' 

\ 

Figure 4.3. Langmuir adsorption isotherm for Model 2 (top) and 
3 (bottom). 



In Figure 4.3 we similarly compare the theoretical isotherm (4.14) with numerical 
values for Model 2 and 3. No unstable cases were detected and the results are all 
in very good agreement. For Model 2, the initial data used in the experiments were 
a bit off the actual equilibrium causing the measured values to "creep" slightly. 
However, all values stay in close proximity to the predicted isotherm curve. 

Finally, in Figure 4.4 the initial approximate surfactant concentration profile is 
compared to the final equilibrium profile for a single value of tp c and for different 
surfactant bulk concentrations ipb- 

Conclusion. We verified the theoretical analysis in Section 3.1 and in particular 
the sharpness of the sufficient condition for ill-posedness (3.5). The Langmuir 
isotherm (4.14) for Model as published originally in [44] and again, using altered 
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-0.5 0.5 -0.5 0.5 -0.5 0.5 

(a) Model 1. (b) Model 2. (c) Model 3. 



Figure 4.4. Initial surfactant concentration profiles according to 
(4.12) (dashed) and numerical equilibrium values (solid) for differ- 
ent values tp b G {1CT 4 , 1(T 3 , 1CT 2 } and ip c = 0.016. For Model 1 
the value of B — B b in (4.13) formally belonging to Model was 
used in order to obtain a closed expression for the initial profile. 

parameters, in [36] could therefore not be obtained. Model 1, although stable in 
most cases, does not a priori possess a natural adsorption isotherm. Model 2 and 
3 are not only stable, but they also satisfy the Langmuir isotherm very accurately. 
Model 3 produces the most sharp surfactant profile and Model 1 yields the most 
damped profile, at least towards the higher values of surfactant concentrations 
tested here. The analytical prediction (4.12) is very accurate for Model 2 and 3. 

4.3. Adsorption dynamics. The adsorption dynamics for an interface in contact 
with a semi-infinite bulk was considered in the early paper [53] by Ward and Tordai. 
Here the time-dependent decrease of surface tension under the presence of a solute 
was explained on the basis of diffusion. Their approach has later been refined; 
notably by Hansen [23] who included also the process of evaporation of solute from 
the interface, and by Ariel, Diamant, and Andclman [9, 10] who proposed a separate 
treatment for ionic surfactants. 

In order to cohere with this classical set-up we continue to use our one-dimensional 
setting with an interface at x = 0. We treat the phase-field variable <f> as being fixed 
in the arguments below; the precise form is not critical as long as the scaling is such 
that the interface is located within \x\ < Cn and such that \<j>(x)\ ~ 1 for |x| 3> Cn. 
As before we assume the constant scalar bulk concentration ipf, = ip(t, x — > oo) and 
further put i^o(t) := ip(t, x = 0), and also V'o.eq := ipo(t °°)- 

The model proposed in [10] is "sharp" in the sense that it is directly formulated 
as a discrete model with a characteristic length-scale a (compartment size) on the 
order of a single molecule. We may therefore refer to such a model as "microscopic" 
and it is clear that any continuous model will have difficulties as the sharpness of 
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the interface approaches this length-scale. Indeed, for the diffuse interface model 
under current consideration, the strict limit Cn — > does not make sense as the 
limiting equilibrium surfactant concentration profile necessarily becomes a constant 
bulk-value with a single point removed at the interface itself (this follows from the 
isotherm relation (4.12)). However, the limiting behavior as Cn — > a <C 1 is still of 
interest since it makes a dynamic comparison with the microscopic model in [10] 
possible. Specifically, we are interested in testing the adsorption dynamics at the 
interface for which there are known discrete modeling approximations available. 

In order to estimate the behavior as Cn — > a, where we may think of a > as 
the length-scale of a single molecule, we write the surfactant concentration in terms 
of inner and outer variables, 

(4.17) ip(x) = ip in (y) for -l<y< 1, with y := x/Cn, 

(4.18) VOO=VW(*) for|a;|>Cn. 

Using the constant phase-field profile <fi(x) — tanh(cc/Cn) we can write the PDE 
(2.2) for V as 

in terms of which 



4 T 4Ex 

(4.20) <p(x) = I (W 2 ) 



(Model and 3) 



Cn V" - f i> + 4 (Model 1 



4 1 4 Ex 



(Model 2) 



For | a; | ^> Cn and in all four cases we readily get the bulk diffusion equation, 

dip out Pi d 2 ip out 



(4.21) 



dt Pe^ dx 2 



Also, through the change of variables y = xj Cn expressed in (4.17), we get the 
inner dynamics 

at Cn Pe^ ay z Cn Pe^, ay 
where for Model 1 the change of variables implies 

(4-23) <p(y) = - 1 J ' V - o ^ + 



2 2 r 4 Ex 

Assuming an asymptotic expansion of the form ip in (t,y) ~ 12j>o Cn J '0i,Cn(^, y) 
we have to leading order in Cn that the inner variable is in equilibrium with the 
outer variable via the Dirichlet (matching) boundary condition ip in (±l) = ipi(t) :— 
ip out (± Cn). As in [9, 10] we may interpret the new variable ip\ as the concentration 
in the sub-surface layer and we have explicitly enforced an even symmetry ip(t, x) — 
ip(t, —x) in order to remain compatible with the modeling in those references. 



18 



S. ENGBLOM, M. DO-QUANG, G. AMBERG, AND A-K. TORNBERG 



We connect the two variables by requiring that the total volume is conserved 
(assuming conservative outer boundary conditions), 

(4.24) 

dV f d^ out f 1 di> ln 2Pi d f 1 

dt J\ x \>cn ot J_ x dt Pe^, dt }_ x 

For Model 2 and 3, with ip in in equilibrium the inner variable can in principle be 
determined explicitly as a function of ip% thanks to an isotherm-like relation. Model 
1 does not possess such an isotherm relation, but in this case a steady-state solution 
of the inner dynamics (4.22) plays the same role. Although an analytic expression 
for the resulting integral is lacking we argue that under the present scaling, the near 
surface dynamics is essentially of discrete character in the two variables ipo and ipi ■ 
As an example, using the trapezoidal rule we get the approximative relation 

Interestingly, (4.21) and (4.25) above correspond directly to Eqs. (2.11) and 
(2.12) in the derivation in [9]. These equations can be solved analytically in terms of 
an integral relation such that asymptotic estimates can be obtained. More precisely, 
when starting with a uniform profile ip(x) — ipb at time t = 0we have the asymptotic 
"footprint" of diffusion [10], 

(4.26) 



There is also an initial transient phase with a \/i-dependence [23] , 
(4.27) ip Q (t) ~ const. + 



n 

Finally, as pointed out in [10], with a uniform initial profile there is also an ultra- 
short linear transient before (4.27) is valid, 

(4.28) - 1 + const, x t. 

4.3.1. Numerical dynamics. This quite interesting typical adsorption behavior is 
displayed in Figure 4.5, again using the numerical method outlined in Appendix B. 
The initial very short linear transient is the phase where the region just outside the 
interface quickly becomes depleted of surfactant. In the following V^-dependcnt 
phase, surfactant diffuses from outside the sub-surface region, enters it, and then 
immediately adsorbs to the interface, keeping the concentration in the sub-surface 
region in equilibrium. The asymptotic regime occurs towards the end of this process 
when the interface gradually becomes saturated and the adsorption therefore slows 
down. 

In Figure 4.6 and 4.7 we numerically fit (using standard polynomial least squares) 
the measured time-dependent interfacial adsorption to the theoretical expressions 
for a sharp interface expressed in (4.26)-(4.28). Although the two models display 
some differences, both versions clearly fit quite well with the theoretical predictions. 
Qualitatively similar results were produced also for Model 1, but the temporal 
scaling is slightly different here. Also, in one of the cases tested, we were unable to 
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Figure 4.5. Adsorption phases. Top: starting from a uniform 
concentration in (a), the interface very quickly drains the sub- 
surface region in (b)-(c). Middle: surfactant now diffuses inwards 
from the bulk while the concentration at the interface steadily 
builds up in (c)-(d). Bottom: in (e), the interface starts to reach 
saturation and in (f), finally, boundary effects due to the finite 
size of the numerical domain become prominent. In this example 
Model 3 was used with parameters Cn = 1/20, ip c = 0.016, Ex = 1, 
and ipb = 10 -2 . Note the logarithmic vertical scale. 



obtain a solution due to instabilities and so we choose not to report these results 
here. 

Conclusion. The Ward-Tordai governing equations (4.21) and (4.25) are approxi- 
mately retrieved by performing a multiscale analysis in terms of inner and outer 
variables. Experimentally, for Model 2 and 3 we verified good agreement as to the 
theoretical behavior in the sharp interface limit. 

5. Numerical method and experiments in 2D 

In this section we perform a final qualitative computational experiment with 
the full hydrodynamic set of equations (2.1)-(2.4). The purpose is to show that 
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1.75 - 




time 

Figure 4.6. Top: initial very short linear transient, and bottom: 
the characteristic diffusion-controlled i^-transient together with 
the £ -1 / 2 -asymptotics. Dashed: numerical values, solid: least- 
squares fit to the theoretical models (see text). These results are 
for Model 2 with parameters Cn = 1/20, ip c = 0.016, Ex = 1. 
Counting from below in the bottom graph we have ipb — [lj 2, 4, 8] x 
10~ 2 (in the top graph the order is the reversed). 



the improved model can capture the nontrivial coupling between the fluid flow 
and the surface forces. This is an essential property as it involves the convective 
redistribution of surfactants, the resulting gradients of surface tension that will 
appear, and a subsequent modification of the overall dynamics. One very generic 
example which we have chosen to study is to show that the presence of surfactant 
in a high enough concentration may inhibit droplet coalescence. 

We present the results from simulations in two spatial dimensions: the extension 
to 3D would involve straightforward extensions of the operators in the model equa- 
tions. With the obvious exception of a computationally heavier solution procedure, 
we do not foresee any other particular differences, since the main physical effects 
of surface stretching, convection of surfactant, diffusion, adsorption and desorption 
are present also in the 2D test cases we have made. 
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2 




time 

Figure 4.7. As in Figure 4.6 but for Model 3 instead. Evidently, 
this model enjoys a slightly faster initial adsorption at the same 
Langmuir-controlled final saturation. 



Given the previous analytical and computational results we restrict the experi- 
ment to Model 3 with chemical potentials given by (3.15)-(3.16). Although we did 
observe the ill-posedness for Model here as well, those results are not detailed 
here. 

5.1. Adaptive finite element method. Successful finite element methods for 
the Cahn-Hilliard equations have been devised by several authors. For a fully 
discrete convergent method targeting the fully coupled hydrodynamic flow, see [29] . 
Versions with degenerate mobilities are more difficult to analyze such that [2, 6] 
are two notable exceptions. Here, although only the no-flow case is treated, on the 
other hand logarithmic terms are allowed in the free energy. 

We carried out our numerical simulations using femLego [12], a software to solve 
general partial differential equations with adaptive finite element methods. The 
PDEs, the boundary conditions, the initial data, and the method of solving each 
equation are all specified in a Maple worksheet such that exact integration and 
differentiation arc possible. 
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The Cahn-Hilliard-type equations (2.1)-(2.2) with chemical potentials according 
to (3.15)-(3.16) are treated as a coupled system for the potentials ^ and fj,^ and 
the composition variables <p an d i 3 - All equations are discretized in space with 
piecewise linear functions and in time using the trapezoidal rule, but with u frozen 
at the previous time-step, thus resembling the strategy in [29]. The coupled system 
of equations is solved using Newton iterations with UMFPACK [8] as the inner 
linear solver. 

To ensure mesh resolution in the vicinity of the interface, an adaptively refined 
and derefined mesh is used with an ad hoc error criterion for each element f2fc, 



The mesh adaptivity is implemented by marking element f2j. for refinement if the 
element size is still larger than the minimum mesh size allowed, h > h m i n , and 
it does not meet the criterion (5.1). In the case that an element does meet the 
criterion, it is marked for derefinement unless it is an element of the initial mesh. 
In the experiment reported here we adjusted TOL and h m - m so that > 10 triangles 
were used across the interfaces, defined here to be the regions where <p varies between 
—0.98 to 0.98. Some more details about this scheme can be found in [12]. 

The Navier-Stokes equations (2.3)-(2.4) are solved using the projection method 
devised in [18]. Again, piecewise linear basis functions are used to discretize space. 
Firstly, an approximate pressure is extrapolated from the previous time-step and 
used when solving the momentum equation. Since this by far is the most costly 
step, we employed an iterative linear solver for this part [17, GMRES, Chap. 3.2]. 
Secondly, a projection step in the form of a Poisson equation derived from (2.3) is 
solved to correct the velocity. 

5.2. Colliding droplets. To demonstrate that our extended model of surfactant 
adsorption can be coupled to hydrodynamics, we have used Model 3 to perform 
simulations of two droplets colliding in a linear shear flow channel. Figure 5.1 shows 
the evolution of the two droplets with droplet Reynolds number Re = 0.5, capillary 
number Ca = 0.1, and surfactant bulk concentration t/>j G {10~ 3 , 10 -4 , 10 -5 }. The 
parameters of the surfactant model for this simulation are given by Cn = 1/20, 
Ex = 0.117, and Pi = 0.5857718. All three cases are plotted on top of each 
other; since increasing the amount of surfactant effectively hinders coalescence of 
the droplets they are easily sorted out. For example, with the smallest value of bulk 
concentration ifib = 10~ 5 , the two droplets have coalesced already at t = 165. For 
V->6 = 10 -4 , they coalesce at t — 190 and they stay separate indefinitely whenever 
ijjf, > 10~ 3 . Due to a comparably strong surface tension effect and a weak fluid 
flow, we do not see a noticeable effect on the droplet's behavior before they are in 
close proximity to each other. However, when two droplets do approach, the fluid 
flow becomes stronger and much more complex as shown in Figure 5.2. 

Figure 5.3 shows a contour plot of the surfactant concentration of two colliding 
droplets in the simple shear flow with ipi, — 10~ 3 . Initially, the increased pressure in 
the gap between the two droplets pushes surfactant away from the near-contact re- 
gion, thus generating a Marangoni stress that affects the droplet-droplet interaction. 
In addition, the reduction of interfacial tension due to the presence of surfactant 
has an effect on droplet deformation. In this way it affects the droplet-droplet 
interaction also as a secondary effect. 



(5.1) 
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Figure 5.1. Evolution of two droplets in linear shear flow in a 
channel. Three cases with surfactant bulk concentration 6 
{1CP 3 , 10~ 4 , 10~ 5 } are displayed here. Top: initial profiles; middle: 
the solutions at time t = 165; bottom: the solutions at t = 190. 
See text for further details. 



6. Summary and conclusions 

In [44], a diffuse interface model for the surfactant adsorption onto the interface of 
two immiscible fluids was presented. Recently, this model was extended to account 
also for different solubility in the two phases, and the possibility to consider systems 
better described by a Frumkin adsorption in addition to the Langmuir adsorption 
[36]. However, the basic model was the same. 
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Figure 5.2. Complex flow (streamlines) for the two approaching 
droplets at t = 75 (ip b = 1CT 3 ). 




FIGURE 5.3. Surfactant concentration on the interface of the two 
droplets in linear shear flow. Left: t — 100; middle: t — 125; right: 
t = 150 (ip b = 10~ 3 ). 



We have shown that this model is mathematically ill-posed for a large set of phys- 
ically relevant parameters, in particular for cases where the interface becomes satu- 
rated with surfactant. We have derived an explicit condition for this ill-posedness, 
and through careful and very accurate one-dimensional simulations illustrated that 
this condition is quite sharp as indicated by instability and blow-up of numerical 
solutions. These conclusions limit the usability of this model. 

We have suggested and analyzed three alternatives to the basic model (Model 
0), as denoted Model 1-3. In Model 1, a natural idea in the form of an energy 
contribution from the gradient of the surfactant concentration is included. However, 
besides from being more complicated, this model cannot reproduce the Langmuir 
adsorption isotherm, and is also found to be numerically unstable in a few cases. In 
Models 2-3, the surface energy accounting for the adsorption of surfactant to the 
interface, is changed from a form containing the gradient of the phase-field function, 
to a "gradient-free form" . This removes the problem of ill-posedness completely, and 
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both models are able to accurately reproduce the Langmuir adsorption isotherm. 
The exact form of this surface energy differs between Models 2 and 3, rendering a 
somewhat different profile for the surfactant concentration. 

The conclusions concerning the mathematical models are also verified in sim- 
ulations with fluid flow. Interaction of droplets in two dimensional shear flow in 
presence of surfactants is taken as a test case. We find that the criteria for well- 
posedness determine the usability of the model also in these more complex cases. 

The phase-field modeling for surfactant laden flows is still in its infancy. In 
this paper, we have highlighted the fact that mathematical well-posedness of the 
resulting equations is not an obvious feature, and that this is an issue that needs 
to be considered as new models are developed. 

6.1. Reproducibility. Our ID spectral-Galerkin code as described in Appendix B 
is available for download at the corresponding author's web-page 1 . Along with it, 
scripts that repeat the numerical experiments in Section 4 are distributed. 
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Appendix A. Non-dimensionalization 

We discuss here in some detail the non-dimensionalization of the diffuse phase- 
field model incorporating surfactants as originally presented in [44]. We follow the 
notation therein closely in what follows. 

Using the same degrees of freedom {(j>, %/}, u} as in Section 2 and with the same 
meaning, but now with <p g [— <pQ, </>q] we have (compare (2.1)-(2.4)) 

0. 

-VP + V ■ (pi/[Vu + (Vu) T ]) - (fN^ - ipVm. 



(A.l) _^ + V -(0u) = 

(A.2) ^+V.(Vm) = 

(A.3) V ■ u = 

(A.4) P (- + u-Vu) = 



Compared to [44], in (A.3)-(A.4) we have changed the setting to incompressible 
flows for which P enforces the incompressibility condition and the two final terms 
in (A.4) acts as surface tension forces [27]. The various free parameters of the model 
are defined as follows. 

Firstly, the free energy still takes the form (2.6) but now in terms of 

(A.5) F , = _^ + | 4 + « (w)2; 

(A.6) F 4 , = kT [i/> log i/> + (1 - V) log(l - V)] , 

(A.7) ^ 1 = _£^(V0) 2 , 

W 

(A.8) F cx = y V0 2 - 



In the scaling of (A.5), the equilibrium planar interface at the origin is given 
by <j>(x) = o tanh(a;/C) with O = y/A/B, ( = y/2n/A, and with {A, B,k} the 
parameters of the Cahn-Hilliard model [27]. 
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The chemical potentials are obtained through variational derivatives of the free 
energy, 

SF 

(A.9) iia, = — = -AS + Bcj) 3 - nAd> + e^Adj + eV-i/> ■ V</> + Wip6, 

ocp 

(A.10) ^ = ^ = ^o gT ^-|(V^ + ^. 

To non-dimensionalize the model, let {</>, ip, u}(i, x) denote the variables and 
coordinates as defined in the preceding paragraph. We scale time and space in 
such a way that (i, x) — (tuo/L,x/L) and change variables according to 4>{t,x) — 
d)o4>(t,x), ij){t,x) = ij){t,x), u(t,x) — UQu(t,x), and also set p — pop, v = v$v. It 
becomes natural to define the new parameters 

(A.ll) Cn = C/i, Ex = e/{W( 2 ), Pi = kT/(A(f) 2 ). 

Note that, as in [44], we tacitly assume that e = k in order to reduce unnecessary 
free parameters. 

The new Ginzburg-Landau energy now becomes (2.6)-(2.10), yielding the chem- 
ical potentials in (2.11)-(2.12). The new definition of the mobilities is 

(A. 12) Afy = M$A, mi/, = m$A<j>l, 

where we recall that a degenerate mobility — 771^-0(1 — ip) is used for ij) in 
(A. 2) (conveniently, — — ip) in (2.2)). It follows that the Peclet numbers 
in (2.1)-(2.2) are obtained as 

(A. 13) Pe^ = Luo/M,f„ Pe^, = Lu /m^,, 

and the Reynolds and Capillary numbers as usually by 

(A.14) Re^, Ca=^°. 

v 7 

Using the closing relation 7 oc k0q/C for the interfacial tension 7 finally yields (2.4). 

Appendix B. Numerical method in ID 

In this section we describe our high-resolution one-dimensional spectral scheme 
for the no-flow case (u = in (2. 1)— (2.2)) , using the original energy (2.7)-(2.10). 
The necessary modifications for Model 1-3 are trivial and are omitted for brevity. 

A brief review of numerical methods for gradient flows in general and surfactant 
models in particular might now be in order. Although many papers treat numerical 
methods for the Cahn-Hilliard and the closely related Allen-Cahn equations, it goes 
without saying that references considering computational surfactant flows are much 
more scarce. In this respect [5] stands out where a hybrid method for surfactant 
two-phase flow is designed. In the setting of a sharp interface treatment and using 
singular perturbation analysis, an integral formulation is used in a semi-implicit 
time-discretization strategy. 

Semi-implicit, or split-step methods otherwise have a long history for gradient 
flows, but is also an active area of research. Tracing their origins back to general 
results in [15, (unpublished)] and [16], two more recent references include [4] and 
[58]. In the former, general fourth order gradient flows are considered, and in the 
latter the Allen-Cahn and Cahn-Hilliard equations are specifically targeted. 
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A somewhat different design strategy is used in [43, 56] where the associated 
Ginzburg-Landau functional is modified in order to provide for numerical stability. 
The design therein relies on (spectral) Galerkin formulations whereas in [4, 37] finite 
differences are rather considered. A finite volume scheme for the Cahn-Hilliard 
equation is obtained in [7] , where the efficiency benefit with time-step adaptivity is 
also stressed. 

On balance, and driven by a need for robustness and transparency rather than 
for efficiency, we have chosen to design a polynomial spectral Galerkin method 
with a very simple but carefully controlled time-discretization strategy. We thus 
postpone more advanced spatial and temporal simulation techniques for another 
occasion. As a word in favor of this set-up, this discretization is similar to other 
Galerkin-based methods such that generalizations to more realistic situations are 
possible. 

B.l. Legendre spectral Galerkin method. We consider semi-discrete weak so- 
lutions (j), ip to (2.1)-(2.2) (u = 0) in the space H 1 ^) = {w; ||w|| + j|Vw|| < oo} 
where O = [-1, 1] and the L 2 (Vt) -inner product and induced norm are understood. 
Using natural (volume preserving) homogeneous Neumann boundary conditions 
and the auxiliary variable $ denoting the chemical potential /lu we get after inte- 
gration by parts the variational formulation 

(B.l) ( x , $) = (x, -4> + tf + ^vV) + ~(Vx, (1 - VW), 

(B.2) ( X; t ) = _^(v x ,V$) 

for all x € H 1 ^) and t € (0,T] assuming available initial data at t = 0. 

The variational formulation for ip is obtained in similar fashion but using explic- 
itly the observation in (2.13). Define first 

Cn 2 - 1 
(B.3) * = -^— (V0) 



4 Ex 



We tacitly assume the boundary conditions Pi dip/dn + ip(l — ip)d^> jdn = on dQ 
and get 

(B.4) ( x , Vt) = - ^-(V X ,Pi V0 + V(l - ^)W), 

Pe^ 

again for all x S i? 1 (il) and t 6 (0, T). In one dimension these boundary conditions 
simplify to 4>"' = (j)' = i/j' = at the endpoints x = ±1. 

As a discrete subspace of H 1 we consider the space of polynomials of degree 
< N and use as test- and trial functions the Legendre polynomials which are or- 
thogonal in the L 2 {[— 1, l])-inner product [40, Chap. 18]. Expanding (f>iq{t,x) — 
T^=o4>n{t)Pn{x) and similarly for ip^ and the auxiliary variable $jv we get from 
the variational form (B.1)-(B.4) the semi-discrete set of equations 

(B.5) M$(t) = a, 

(B.6) M4>'(t) = a, 

(B.7) My/(i) = /3. 
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The mass-matrix M is diagonal with entries Ma = \\Pi\\ 2 = 2/(2i + 1) for i E 
{0, . . . , N} and 

(1 \ Cn^ 

Pi, (pN + 4?N + ^^n4>n \ + ~^( P i, (! ~ *I>nWn), 

(B.9) ai = -pL(P/,^), 

(B.10) A = -~(P/,Pi^ + ^Jv(l - 

with ^ defined by (B.3) (hence imposed strongly). The derivatives are readily 
obtained by explicit differentiation of the Legendre expansions [25, Appendix B.1.4]. 

To evaluate the inner products we use the associated Gauss-Legendre quadra- 
ture of at least the same order as the scheme, often considerably higher in order to 
avoid or at least mitigate aliasing errors [25, Chap. 6]. It is known that, for non- 
linear operators, aliasing errors and Gibbs oscillations near sharp gradients may 
drive the scheme unstable even for smooth problems [24]. This problem is perhaps 
mainly associated with shocks but is of relevance also here since phase-field models 
naturally produce sharp gradients in the vicinity of the interface. 

We have therefore implemented filtering along the lines presented in [25, Chap. 9.2]. 
To this end we replace the Legendre expansion with the filtered expansion 

N 

(B.ll) <MM) = ^^(n/N)$ n (t)P n (x), 

n=0 

and similarly for ipjy and $^r. A suitable filter is the exponential one defined by 
[25, p. 164] 



(B.12) 0(r,) 



1 V < Vc 



with [rj c ,p, a] filter parameters. Following the recommendations in [24] we take the 
order to be p — 12 and set a = — log e « 36 with e the double precision machine 
accuracy. We also specify r/ c = 0.25 so that the 25% low modes are not filtered 
at all. The filtered method is then conveniently implemented by replacing the 
mass-matrix in (B.5)-(B.7) with M defined by Ma := ^{i/N^Mu. 

B.2. Adaptive discretization of time. Since we are considering problems that 
either develop instabilities or for which there are nearby ill-posed problems, rather 
than aiming for efficiency and high order we simply use the backward Euler method 
for discretizing (B.5)-(B.7) in time: 

(B.13) M$ fe+1 = a k+1 , 

(B.14) M((f> k+1 -4> k ) = Ata k+ \ 

(B.15) M{4> k+l - V?) = At/3 fe+1 , 

with At a suitable time-step and where the dependence on time is expressed through 
superscripts. 

Although there are arguably more efficient split-step methods available that 
would probably apply here [4, 37, 58], the backward Euler method has a specific 
advantage for gradient flows such as the one at hand. Namely, it is one of very few 
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unconditionally gradient stable methods (if not the only one), which preserves decay 
in the exact energy functional of the system (cf. [49, Chap. 5.6.1]; Theorems 5.6.1- 
5.6.3). Correct energy dissipation has been recognized as one of the most important 
criteria when designing numerical methods targeting gradient flows, which in spe- 
cific cases may display chaotic behavior and "weak turbulence" (see [48, Sect. 4] 
and [7]). An intuitive argument is the fact that the Ginzburg-Landau energy is 
the modeling step and controlling the rate of energy dissipation is therefore more 
important than formal order. 

For the Cahn-Hilliard equation it has been demonstrated by example that time- 
step adaptivity is of vital importance since the dynamics is rich in scales and since 
rapid transients due to initial data and boundary conditions occur naturally [7, 37] . 
We address this by using error estimates from half-steps and a careful time-step 
selection mechanism. 

Modern adaptivity based on digital filtering was used in our code [45, 46]. To 
discuss it and for brevity, let Ta< be the forward-in-time map by the backward 
Euler method as obtained by solving (B.13)-(B.15). We put u\ = T^tUo and 
U2 = ?At/2 u o f° r some initial state uq. By Richardson extrapolation this gives us 
a point-wise error estimate 

(B.16) err = -1*21/3. 

For given relative and absolute tolerances Rtol and Atol we compare the error to the 
effective tolerance vector TOL := max(Rtol |u 2 |, Atol) through the control variable 

(B.17) cerr := L(|| err / TOL ||^ 1/2 ), 

(element-wise division), where 1/2 is the order of the error estimate and where L is 
a denoising limiter (see below) . We consistently used the fairly stringent tolerances 
Rtol = 10~ 6 x diag(M -1 / 2 ) (vector) and Atol = 10~ 8 in our experiments. 

The control objective is now to keep cerr = 1. An updated step-size is obtained 
from a previous time-step At' via the relation At — pAt' , where in terms of 

(B.18) p := L(F(cerr, cerr', p')) 

and where primes denote variables computed at the previous time-step. As to the 
choice of F, the step-size controller, we decided upon H211b defined by 

(B.19) F(a,P,j) =a 1/b p 1/b j 1/b 

with the recommended choice of the single parameter being 6=1/4 [45]. Both 
(B.17) and (B.18) employ the arctan- limiter [46, Sect. 6] 

(B.20) L(a) = 1 + k arctan 

with, respectively, n = 2 in (B.17) and n = 1 in (B.18). The limiter's purpose in 
the latter case is to smoothly restrict the rate by which the step-size may change 
and, in the former case, to reduce the impact of noise in the error estimate. 

By the very action of the step-size controller, the sequence of ps is generally 
much smoother than the corresponding values for the control variable in (B.17). 
For computational stability and following a suggestion in [46], we therefore base 
the decision of rejecting steps on the former variable rather than on the latter. 
Specifically, given some tolerance factor M (we took M = 2), the controller is 
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scheduled to reject the step whenever err > M ■ TOL. Following the steps that lead 
to (B.18) this suggests the rejection criterion 

(B.21) p < Pmin := L(F(L(M~ 1 ^ 2 ), 1, 1)) « 0.9179. 

In practice and with this strategy for time-step selection, step rejections occur very 
rarely. 

For solving the nonlinear system of equations (B.13)-(B.15), simplified Newton 
iterations with a numerical Jacobian were used. We followed the prescription for 
estimating the speed of convergence in [21, Chap. IV. 8], but also found it rewarding 
to augment the code with control strategies devised in [20] for reusing previous 
Jacobian evaluations and factorizations. 

The above described code was used successfully over a wide range of test-cases. 
For reproducibility, the Matlab source code is made freely available from the corre- 
sponding author's web-page. In particular, scripts generating all numerical results 
in Section 4 are distributed (see Section 6.1 for further information). 
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